function seiscell = fmodelFastAVA(type, DenModel, VpModel, VsModel, angles, wavelet)
%% FORWARD MODEL FOR REALIZATIONS 
%
% [input]
%   DenModel [grid]      - DENSITY REALIZATION IN GRID
%   VpModel  [grid]      - VP REALIZATION IN GRID
%   VsModel  [grid]      - VS REALIZATION IN GRID
%   angles [1 x angles]  - STACKS OF PARTIAL ANGLES 
%   wavelet              - WAVELET GATHER USED FOR CONVOLUTION
%   nullDataValue        - NULL VALUES FROM ORIGINAL DATA THAT WERE MASKED 
%
% [output]
%   seiscell             - CELL OF COMPUTED SYNTHETIC SEISMIC
%
% [functions]
%   convn.m
%
% [MODIFICATIONS]
%   JUL 2019 - update on code match slow version. handles NaN
%   JAN 2017 - change to match BF convention.
%
% LA - JAN 2017
%

%% ALLOCATE AND INTIALIZE
seisgrid    = zeros(size(DenModel,1) * size(DenModel,3), size(DenModel,2) * size(angles,2));
switch type
    case 1 % Aki-Richards
        % COMPUTE VARIABLES FROM DENSITY, Vp AND Vs MODELS
        Dden = diff(DenModel,1,3);
        DVp  = diff(VpModel,1,3);
        DVs  = diff(VsModel,1,3);

        den  = (DenModel(:,:,1:end-1) + DenModel(:,:,2:end))./2;
        Vp   = (VpModel(:,:,1:end-1)  + VpModel(:,:,2:end)) ./2;
        Vs   = (VsModel(:,:,1:end-1)  + VsModel(:,:,2:end)) ./2;
        for n = 1:size(angles,2)

            % angles in degrees
            ct   = cosd(angles(1,n));
            p    = sind(angles(1,n)) ./ Vp;

            Rpp                 =  zeros(size(DenModel,1), size(DenModel,2), size(DenModel,3));
            Rpp(:,:,1:end-1)    = (0.5 .* (1 - 4 .* p.^2 .* Vs.^2) .* (Dden ./ den)) + ( DVp./(2 .* ct.^2 .* Vp  )) - ( 4 .* p .^2 .* Vs .* DVs); 

            % convert NaN in Rpp to zero
            Rpp(find(isnan(DenModel) == 1)) = 0;

            % convolve
            temp                = permute(convn(permute(Rpp,[3 2 1]),wavelet(:,n),'same'),[3 2 1]);

            % set back NaN
            temp(find(isnan(DenModel) == 1)) = NaN;

            seisgrid(:,((n-1) * size(DenModel,2) + 1): ((n-1) * size(DenModel,2)) + size(DenModel,2)) = reshape(permute(temp,[3 1 2]), [], size(DenModel,2));
        end
    case 2 % full zoeppritz       
        for n = 1:size(angles,2)
            Rpp                 =  zeros(size(DenModel,1), size(DenModel,2), size(DenModel,3));
            angleRadians        = angles(1,n) * (pi./180.0);  
            p                   = sin(angleRadians) ./ VpModel(:,:,1:end-1);
            
            Dden = DenModel(:,:,1:end-1);
            DVp  = VpModel(:,:,1:end-1);
            DVs  = VsModel(:,:,1:end-1);
            DVsSq    = DVs .^ 2;
            DVsP     = DVs.*p;
  
            den  = DenModel(:,:,2:end);

            Vp   = VpModel(:,:,2:end);
            VpSq    = Vp .^ 2;
            
            Vs      = VsModel(:,:,2:end);
            VsSq    = Vs .^ 2;
            VsP     = Vs.*p;
            
            Rpp(:,:,1:end-1)    = ((((den .* (1. - 2.*(VsP).^2)) +...
                2 .* Dden .* (DVsP).^2) .* ((cos(angleRadians)) ...
                ./ DVp) - ((Dden .* (1. - 2.*(DVsP).^2))...
                + 2 .* den .* (VsP).^2) .* (sqrt(1./(VpSq) - p.^2)))...
                .* (((den .* (1. - 2.*(VsP).^2)) + 2 .* Dden...
                .* (DVsP).^2) .* (sqrt(1./(DVsSq) - p.^2)) + ((Dden...
                .* (1. - 2.*(DVsP).^2)) + 2 .* den .* (VsP).^2)...
                .* (sqrt(1./(VsSq) - p .^2))) - (((den .* (1. - 2.*(VsP).^2))...
                - (Dden .* (1. - 2.*(DVsP).^2))) + (2 .*(den .* Vs...
                .^2 - Dden .* DVs .^2)) .* ((cos(angleRadians)) ./ DVp)...
                .* (sqrt(1./(VsSq) - p .^2))).* (((den .* (1. - 2.*(VsP).^2)) ...
                - (Dden .* (1. - 2.*(DVsP).^2))) - (2 .*(den .* Vs...
                .^2 - Dden .* DVs .^2)) .* (sqrt(1./(VpSq) - p.^2))...
                .* (sqrt(1./(DVsSq) - p.^2))) .*p.*p)./((((den...
                .* (1. - 2.*(VsP).^2)) + 2 .* Dden .* (DVsP).^2)...
                .* ((cos(angleRadians)) ./ DVp) + ((Dden .* (1. - 2.*(DVsP).^2))...
                + 2 .* den .* (VsP).^2) .* (sqrt(1./(VpSq) - p.^2))) .* (((den...
                .* (1. - 2.*(VsP).^2)) + 2 .* Dden .* (DVsP).^2)...
                .* (sqrt(1./(DVsSq) - p.^2)) + ((Dden .* (1. - 2.*(DVsP).^2))...
                + 2 .* den .* (VsP).^2) .* (sqrt(1./(VsSq) - p .^2)))...
                + (((den .* (1. - 2.*(VsP).^2)) - (Dden .* (1. - 2.*(DVsP).^2)))...
                - (2 .*(den .* Vs .^2 - Dden .* DVs .^2))...
                .* ((cos(angleRadians)) ./ DVp) .* (sqrt(1./(VsSq) - p .^2)))...
                .* (((den .* (1. - 2.*(VsP).^2)) - (Dden .* (1. - 2.*(DVsP).^2)))...
                - (2 .*(den .* Vs .^2 - Dden .* DVs .^2))...
                .* (sqrt(1./(VpSq) - p.^2)) .* (sqrt(1./(DVsSq) - p.^2))) .* p .* p);
        
            Rpp(:,:,end)        = 0;  
            
            Rpp = real(Rpp);
            % convert NaN in Rpp to zero
            Rpp(find(isnan(DenModel) == 1)) = 0;

            % convolve
            temp                = permute(convn(permute(Rpp,[3 2 1]),wavelet(:,n),'same'),[3 2 1]);

            % set back NaN
            temp(find(isnan(DenModel) == 1)) = NaN;

            seisgrid(:,((n-1) * size(DenModel,2) + 1): ((n-1) * size(DenModel,2)) + size(DenModel,2)) = reshape(permute(temp,[3 1 2]), [], size(DenModel,2));
        end        
end
% index to resort seismic in cell
idx = zeros(1,size(seisgrid,2));
for n = 1:size(DenModel,2)
    idx((n-1)*(size(angles,2))+1:n*(size(angles,2))) = n:size(DenModel,2):size(DenModel,2)*size(angles,2);
end
seisgrid = seisgrid(:,idx);

seiscell = mat2cell(seisgrid, repmat(size(DenModel,3),[size(DenModel,1),1]),  repmat(size(angles,2),[size(DenModel,2),1]));

end